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wave  excitations  on  a  large  pontoon  and  tests  the  feasibility  of  reducing  a 
nonlinear  free  surface  problem  by  perturbation  expansions.  A  simulation 
model  has  been  developed  based  on  the  perturbation  expansion  technique 
to  estimate  the  wave  forces.  The  model  uses  a  versatile  finite  element 
procedure  for  the  solution  of  the  reduced  linear  boundary  value  problems. 
This  procedure  achieves  a  fair  compromise  between  computation  costs 
and  physical  details  by  using  a  combination  of  2D  and  3D  elements.  A 
simple  hydraulic  model  test  was  conducted  to  observe  the  wave  forces 
imposed  on  a  rectangle  box  by  Cnoidal  waves  in  shallow  water.  The  test 
measurements  are  consistent  with  the  numerical  predictions  by  the  simula¬ 
tion  model.  This  result  shows  favorable  support  to  the  perturbation  ap¬ 
proach  for  estimating  the  nonlinear  wave  forces  on  shallow  draft  vessels. 
However,  more  sophisticated  model  tests  are  required  for  a  full  justifica¬ 
tion.  Both  theoretical  and  experimental  results  show  profound  second- 
order  forces  that  could  substantially  impact  the  design  of  ocean  facilities. 
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that  the  amplitudes  of  waves  and  structure  motions  are  small  in  comparison  with  the  wavelength 
and  the  nominal  dimension  of  structures.  Thus,  the  diffraction  waves  from  the  structure  may  be 
approximated  by  a  simple  harmonic  function.  These  assumptions  preclude  their  application  in 
severe  sea  states.  Large  waves,  in  fact,  present  a  much  steeper  profile  than  that  of  a  simple 
harmonic  function.  A  precise  description  of  the  wave  kinematics  requires  retention  and  proper 
treatment  of  higher  order  terms  of  the  equation  of  motion  and  a  proper  specification  of  the  water 
body  boundary  at  the  free  surface.  These  two  factors,  which  are  ignored  in  linearized 
approximations,  are  the  principle  sources  of  wave  nonlinearity.  When  these  nonlinear  effects  are 
included  in  the  diffraction  of  waves  by  a  body  there  are,  at  second  order,  interactions  at  the  sums 
and  differences  of  the  component  frequencies  of  the  incident  waves.  Although  the  magnitudes 
of  these  nonlinear  effects  are,  in  general,  only  second  order,  they  act  at  frequencies  away  from 
that  of  the  ambient  wave  energy,  and  may  therefore  be  of  primary  concern.  This  is  especially 
true  when  these  excitations  are  near  the  natural  periods  of  the  body  motions  or  where  restoring 
or  damping  forces  are  small.  With  the  presence  of  large  waves,  second-order  effects  are 
important  corrections  to  the  linearized  results. 

Theoretical  developments  of  the  second-order  diffraction  problem  have  until  recently  been 
limited.  Studies  on  simple  geometry  of  the  uniform  vertical  circular  cylinder  were  conducted 
by  Issacson  (Ref  1),  Chakrabarti  (Ref  2),  Molin  (Ref  3),  Wehausen  (Ref  4),  Hunt  and  Baddour 
(Ref  5),  Chen  and  Hudspeth  (Ref  6),  Rahman  (Ref  7),  and  Ogilvie  (Ref  8).  The  results  are, 
however,  controversial  due  to  the  difficulties  in  the  correct  treatment  of  the  second-order  free 
surface  boundary  condition  and  a  proper  specification  of  the  radiation  condition  for  the  second- 
order  diffracted  waves.  Previous  studies  indicate  that  a  weakly  nonlinear  diffraction  problem 
may  be  reduced  to  a  set  of  linear  boundary  problems  by  the  perturbation  method.  Solutions  may 
be  obtained  by  Green’s  function  method  or  the  finite  element  method.  Kim  and  Yue  (Ref  9) 
presented  a  second-order  diffraction  solution  for  an  axisymmetric  body  in  constant  water  depth 
following  the  former  approach.  The  solution  gives  an  explicit  second-order  flow  potential  in 
terms  of  wave-source  Green  functions.  The  nonlinear  forces,  accurate  to  the  second  order, 
include  slowly  varying  drift  forces  and  biharmonic  forces,  which  appear  at  a  frequency  twice  that 
of  the  incident  wave.  The  present  study  pursues  the  finite  element  method  approach,  which  gives 
the  possibility  of  addressing  a  sloping  sea  bottom  of  irregular  geometry. 

An  immediate  need  for  these  analytical  techniques  exists  at  the  Naval  Civil  Engineering 
Laboratory  (NCEL)  in  the  design  of  connectors  for  a  pontoon-supported  waterfront  facility.  The 
requirement  of  including  a  second-order  diffraction  force  in  the  assessment  of  extreme 
environmental  forces  imposed  by  maximum  design  waves  is  an  apparent  trend  in  the  near  future. 
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PROBLEM  DEFINITION 


A  mathematical  idealization  of  hydrodynamic  couplings  between  waves  and  a  group  of 
large  floating  or  fixed  structures  in  water  with  a  free  surface  is  defined  in  Figure  1.  Wave 
activity  around  the  structures,  hydrodynamic  loads  on  the  structures,  and  the  resulting  motion 
of  the  structure  are  investigated.  The  water  body  may  be  open  or  partially  sheltered,  and  the 
water  depth  may  vary.  Relevant  waves,  including  incident,  scattering,  and  radiation  waves,  are 
described  in  a  way  similar  to  Stokes’  second-order  theory.  The  water  body  is  divided  into  a 
three-dimensional  (3D)  region  D, ,  and  a  two-dimensional  (2D)  region  Q2.  The  3D  region 
includes  areas  near  structures  where  the  seafloor  may  change  substantially  and  the  2D  region 
includes  the  remaining  areas  where  the  seafloor  is  reasonably  flat.  In  3D  regions,  waves  are 
addressed  in  detail  in  3D  space  to  account  for  the  effects  due  to  irregular  geometries  of  structures 
and  the  seafloor.  In  2D  regions,  waves  are  approximated  by  plane  wave  theories  (Ref  10)  in 
terms  of  free  surface  parameters  for  better  computation  efficiency.  Boundaries  of  these  two 
regions  are  illustrated  in  Figure  1 .  These  boundaries  are  the  radiation  boundary,  ?2 ,  which 
truncates  the  water  domain  to  a  finite  area,  and  the  border  separating  the  two  regions,  F, . 
Waves  transmitted  beyond  the  radiation  boundary  are  accounted  for  in  a  collective  form  in  terms 
of  parameters  on  the  radiation  boundary.  Incident  waves  are  specified  on  the  wave  input 
boundary  T, ,  which  coincides  with  the  radiation  boundary  for  convenience,  when  wind  wave 
excitations  are  considered.  Other  boundaries  involved  are  the  sea  bottom,  rb,  the  wetted 
surfaces  of  floating  structures,  T§ ,  and  fixed  structures,  Tw .  Subdivisions  of  the  water  body  and 
the  associated  boundaries  are  illustrated  in  Figure  2. 


PROBLEM  FORMULATION 

A  Cartesian  coordinate  system  (OXYZ)  fixed  to  the  earth,  as  shown  in  Figure  3,  is 
employed  as  an  inertia  reference,  where  the  X-Y  plane  lies  on  the  mean  water  surface  and  the 
Z  coordinate  is  measured  positive  upward  from  the  still  water.  The  responses  of  a  floating  body 
are  described  in  reference  to  a  body  coordinate  system,  Gxyz,  attached  to  the  body,  with  its 
origin  G  located  at  the  center  of  gravity  of  the  body.  The  xy  plane  is  parallel  to  the  static 
waterplane  and  the  z  coordinate  is  directed  vertically  upward.  For  the  purpose  of  computing 
wave  drift  forces  and  moments  exerted  on  the  body,  a  third  system  of  coordinate  axes  (Gxyz) 
is  defined  where  its  origin  is  at  the  center  of  gravity  of  the  body  and  is  always  parallel  to  the 
fixed  coordinate  system  OXYZ. 

With  the  assumption  that  the  fluid  is  inviscid  and  incompressible  and  that  the  flow  is 
irrotational,  the  relevant  wave  problems  can  be  addressed  using  potential  flow  theory. 

Assuming  weakly  nonlinear  waves,  the  total  velocity  potential  *  is  expressed  as  a 
perturbation  in  the  wave  slope  parameter  e  e  kA<  1: 

=  e  +  e2  +  0(e3) 
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where: 

#0)  =  first-order  velocity  potential 
*(2>  =  second-order  velocity  potential 

Further,  k  is  the  incident  wave  number  given  by  the  dispersion  relationship  o2  =  gk  tanh  (kh), 
and  g  is  the  gravitational  acceleration.  The  associated  free  surface  profile  rj  is  decomposed 
likewise: 

q  =  eq(,)  +  +  0(e3) 


This  approximation  leads  to  the  decomposition  of  the  otherwise  nonlinear  boundary  value 
problem  to  a  set  of  two  linear  boundary  value  problems.  The  first  order  problem  has  been 
described  in  Reference  1 1 .  This  discussion  will  be  limited  to  the  second  order. 

The  first-  and  second-order  potentials  must  comply  with  the  Laplace  equation  and 


appropriate  boundary  conditions  as  follows. 

First-Order  Boundary  Value  Problem 

Laplace  equation  (fluid  domain):  V2  *  0 

(2) 

Seabed  and  fixed  structure  boundary:  -  = 

dn 

0 

(3) 

Floating  structure  boundary:  -  =  vn 

on 

(4) 

Free  surface:  g  =  0 

(5) 

Radiation  condition:  lim  y/R  [  — —  -  i  k  ] 

r—  (3R  ) 

=  0 

(6) 

Second-Order  Boundary  Value  Problem 

Laplace  equation  (fluid  domain):  V2  00)  =  0 

(7) 

d 

Seabed  and  fixed  structure  boundary:  -  = 

dn 

0 

(8) 

Floating  structure  boundary: 
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d<X>G) 

For  radiation  waves:  - 

dn 


(9) 


For  scattered  and  incident  waves: 

—  =  -(x(1)-V)V$(1)*n  +  (vO)-v*®)-Nfl)  at  body  surface.  (10) 

dn 

Free  surface: 

+  g  +  on  Z  =  0.  (11) 

Radiation  condition:  outgoing  waves  from  the  floating  bodies. 

In  the  above,  R  =  (x2  +  y2)1'2  is  the  radial  distance  from  the  origin  of  the  inertia  frame, 

and  —  is  the  normal  derivative  into  the  body.  The  second-order  problem  is  complicated  by 
dn 

the  inhomogeneous  forcing  term  in  the  free  surface  boundary  condition  (Equation  10),  which  is 
given  in  terms  of  quadratic  products  of  the  first-order  potential.  The  first-order  problem  is 
classical  and  a  variety  of  numerical  solutions  are  available.  For  example,  a  solution  using  a 
finite  element  scheme  has  been  obtained  in  a  previous  study  (Ref  11).  Equations  10  and  11 
correlate  with  Hence,  is  completely  defined  once  is  known. 

The  first-order  velocity  potential  may  be  expressed  in  component  waves  as  follows.: 

*0)  =  oj1)  +  +  #jj>  (12) 


where: 

=  first  order  incident  wave  potential 
»  fust  order  diffraction  wave  potential 
♦{J*  =  first  order  radiation  wave  potential 


For  monochromatic  incident  waves,  the  time  dependency  is  separated  and  rewritten: 

*<»  .  [<  *  03) 
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Substituting  Equation  13  into  the  second-order  free  surface  boundary  condition  in  Equation  11, 
a  quadratic  periodic  forcing  function  is  obtained  which  oscillates  at  twice  the  frequency  of  the 
linear  waves.  As  a  result,  Equation  1 1  may  be  written  as: 


=  Q(x,y,0)e-2iot  <14) 

where  the  quadratic  forcing  function  Q(x,y,0)  consists  of  nine  components  as  follows: 

Q(x»y»0)  =  Qn  +  Qss  +  Qbb  +  ^is+Qsi^  +  (Qib  +  Qbi^  +  (QsB+Qbs) 

where: 

Qjj  =  plane  wave  forcing  function  which  leads  to  the  ordinary  Stokes’ 

second-order  wave  component  at  double  frequency  2o> 

QsS’Qbb  =  interactions  of  first-order  scattered  waves  and  first-order 
radiation  waves,  respectively 

QjS,...,Qbs  =  cross  interaction  of  the  first-order  terms 
The  details  of  these  components  are  given  in  the  following: 

Qji  -  l  CO  2(V cpj  y  -  <Pj  $izz  +  — 9i  tiz 

g 

Qss  -  i<of2(V^V  -  4>s }  4>szz  ♦  ^4?  *8 

l  g  J 

n  {„  x(1)  x(1)  .  o>2  .(1)  .  (1) 

Qbb  =  lw  2<v<»>b  )  -  9b  $bzz  +  — 9b  9bz 

g 

ft  j  x(l)  j.0)  .  ti)  j,(l)  i(l) 

Q.e  -  iw  2(V9i  )(V9s  )  ~  9j  9szz  +  9i  9sz 

g 
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‘S  I 


=  10) 


2(V*?,)(V*j‘’)  -  4>s j  4>izz  *  —  Oiz 

g 


k«>) 


«  .k<»>  A<» 


(16) 


£1B 


=  1  0) 


2(V4>11))(V441))  -  ♦?>♦&*  ♦  —  4h0> 4>bz 

g 


Qbi 


1  0) 


2(Vtf)(Vtf)  -  ♦S)4>gz 


<t4J>4> 


(1) 

IZ 


CSB 


=  10) 


2(V<d<1>)(V44,))  -  4>S } 4>bzz  - 

g 


CBS 


=  10) 


2(V<|4,,)(V*<I>)  -  44" Viz  ♦  —  C+'s" 

g 


where  the  subscriptions  of  spatial  coordinates  indicate  partial  differential.  These  quadratic 
forcing  function  terms  behave  like  a  nonuniform  pressure  field  applied  to  the  free  surface  in  first- 
order  problems  as  described  by  Wehausen  and  Lai  tone  (Ref  12).  They  generate  additional 
cylindrical  standing  and  outwardly  going  progressive  waves  at  the  double  frequency. 

Since  Equation  14  is  nonhomogeneous,  the  second-order  potential  consists  of  the 
homogeneous  solution  *(2>H  and  the  particular  solution  **2lp  as: 

#<2>  =  $<2)P  +  $<2)H 

(17) 


These  components  satisfy,  respectively,  the  homogeneous  and  inhomogeneous  force  surface 
condition  (Equation  11),  and  jointly,  the  inhomogeneous  body  boundary  condition.  Equation  9 
or  10.  The  particular  solution  may  be  further  separated  as: 

4>a)P  =  *u  ♦  *?s  +  *S»  *  (®?s  +  *?i) 

(18) 

♦  (*2  +  *2) +  (•£ +  *?b) 
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The  homogeneous  solution  consists  of: 


$<2>H  =  <*,£>  + 


By  substituting  Equations  18  and  19  into  Equation  17,  the  total  second-order  potential  becomes: 

-  *ff*[*S**?.*<*?.**S>*<*8**8> 

(20) 

+  (*bs  +  *?b)]  + 


®a)  =  *5 ?  ♦  o£s  + 


where: 


-  [*S  -  *?.  ♦  (*s  ♦  *S)  *  <*2  *  *S)  -  (*s  *  *?,)]  -  •? 


The  component  represents  the  forced  wave  motion  generated  by  the  Qn  forcing  function, 

which  is  the  usual  second-order  Stokes’  wave.  is  the  second-order  scattered  wave  potential 

and  is  the  second  order  radiation  wave  potential  generated  by  body  motion.  Both  satisfy 
the  homogeneous  free  surface  boundary  condition.  Their  far  field  behavior  is  given  by: 


4>H  -  r  ^  t**  +  0  (r  2),  r  >  1 


where  k2  is  the  double  frequency  wave  number  satisfying  4a>2  =  gk2  tanh  (kh).  This  means  that, 
in  principle,  these  components  can  be  solved  by  the  existing  first-order  method. 

The  second-order  potential  *(2)  in  Equation  21  consists  of  three  parts:  (1)  the  incident 

wave  component  G®,  (2)  the  pseudo-diffraction  s,  and  (3)  the  radiation  components 
The  second-order  incident  wave  potential  is  given  by: 
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<x> 


a> 

ii 


2  ft  H  /  rc  H\  cosh[2k(d  +  z)]  e_i2<0t 
8  kT  [,  L  ;  sinh4(kd) 


(22) 


The  radiation  components  satisfy  the  same  equations  as  the  radiation  potentials  of  the 
first-order  potential  theory,  at  the  double  frequency  2o>.  Newton’s  law  applies  through  the  use 
of  Bernoulli’s  equation  in  the  computation  of  the  second-order  wave  exciting  loads,  which  results 
in  the  equations  of  motion  at  the  second  ord't,  to  be  discussed  later. 

The  solution  procedure  for  Equation  14  using  finite  element  method  is  straightforward, 
regardless  of  the  complexity  of  the  forcing  function  on  the  right-hand  side  of  the  equation. 


FINITE  ELEMENT  FORMULATIONS 

According  to  the  calculus  of  variations  (Ref  13),  the  solution  to  a  boundary  value  problem 
is  the  potential  which  minimizes  a  certain  functional,  in  terms  of  the  governing  equation  of  the 
problem.  For  *(2),  the  functional  n(2)  can  be  expressed  as: 

n®  =  [  ( f  1  (V4>(2))2dQ  +  f  f  fOfr^dS  (23) 

J  J  *  2  ^  + 

c 

where: 

n  =  fluid  domain 

Sf  =  free  surface 

SB  =  body  boundary 

SR  =  boundary  for  the  radiation  condition 

The  integral  (//  )  over  the  unknown  free  surface  in  Equation  23  may  be  approximated 

by  a  Taylor  series  expansion,  in  terms  of  the  variables  at  the  mean  sea  surface. 

Solutions  for  the  second-order  diffraction  and  radiation  boundary  value  problems  may  be 
formulated  using  the  standard  finite  element  method.  The  radiation  boundary  condition  is 
imposed  on  a  finite  boundary,  r2,  at  moderate  distances  from  the  structures  of  interest.  The 
entire  fluid  inside  this  boundary  is  divided  into  two  regions.  Each  region  has  the  functionals 
corresponding  to  its  respective  governing  equations  and  associated  boundary  conditions  as 
follows: 
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For  the  three-dimensional  region: 


n(2>  =  Iff  i  (V<tf)2dQ  -  / / 2  —  («*>f >)2 dS 

o,  *,  8 

_  1 1  4f>dS  -  // (24) 

s,  S  s, 

♦  / /  [(x(1)  *V)  V <J>(1)  fi  +  (v(1)  -  V 4><!))  -N(,)]im  •  6tan4>fk)dSk 

s. 


j  =  II,D+S,i,k,... 
i  =  1,2, ...6 
k  =  1,2,.  ..M 

where  M  is  the  total  number  of  floating  bodies  and  4>ik  is  that  of  the  radiated  waves  due  to 
motion  in  the  i0*  mode  by  the  k*  body,  in  which  i=l,2,...,6  corresponds  to  the  surge,  sway, 
heave,  roll,  pitch,  and  sway  modes,  respectively.  The  quantity  Q  in  the  third  term  in  the  above 
equation  is  given  in  Equation  15. 

For  the  two-dimensional  region: 

D®  =  J/-{F*(V<|>f))2-4k2F,(<t>f))2}dQ2 
°2  2 

(25) 

dd>C) 

-  / F*  —-2-  ifrPdS  ♦  (Radiation  Condition  for  j  =II,D+S,i,k... 

i  dn 


The  function  F*  in  Equation  25  is  the  function  defined  by  Equation  A-22  of  Appendix  A. 
On  Tj,  the  three-dimensional  region  is  coupled  with  the  two-dimensional  region  by  the 
following  relationship: 


4>j0)(x,y^)  *  ^(x.y) 


cosh  2  k(d  +  z) 
sinh4(kd) 


j=H,S+D,i,k 


(26) 


where  d  is  the  water  depth. 
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FINITE  ELEMENT  DISCRETIZATION 

The  integral  Equations  24  and  25  define  the  velocity  potentials  of  the  water  flow  over  the 
entire  fluid  domaii .  The  fluid  and  boundaries  are  thus  discretized,  as  shown  in  Figure  2,  to 
conform  the  finite  element  applications. 

The  potential  within  each  element  can  be  approximated  in  terms  of  the  unknown 
nodal  values  { <J>^  }  of  that  element  and  interpolation  functions  [N].  In  matrix  notation: 

4>®'  =  [N]{<|>?}  (27) 


The  functional  is  now  minimized  with  respect  to  the  nodal  values  <J>®  for  a  single 
element  to  determine  the  flow  potential  which  satisfies  the  boundary  value  problem.  For  the  3- 
dimensional  case,  Equation  27  can  be  inserted  into  Equation  24  and  written  as: 


an®  - 


an®  „  .  (2> 


a*; 


(2) 


s  4>r  =  o 


This  being  true  for  any  variation,  6<t>®  requires: 


an® 


a<t>e 


(2) 


fan®/d(<l>?)j 

an®/0(4>?)2 


an®/0(<l>®), 


=  o 


(28) 


(29) 


from  which  parameters  are  to  be  determined.  Because  the  functional  is  quadratic,  element 
Equation  29  reduces  to  a  standard  linear  algebraic  equation,  that  is: 


an® 

a*? 


tki {<*,?}  -  in  =  o 


(30) 


This  transforms  the  integral  equations  into  a  set  of  linear  algebraic  equations  in  the  unknown 
value  of  velocity  potential  at  boundary  nodes  on  all  elements. 
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Stiffness  Matrix  and  Nodal  Force  Vector  for  the  3D  Element 

The  element  stiffness  matrix  [k]  and  the  element  nodal  force  vector  {f}  in  Equation  30 
are  given  as: 

[k]  =  fff  (VN)T(VN)dxdydz  -  f  f  —  NTNdS  (31) 

(O),  <sf>.  8 


|f)  »  -/ /  NT(vft)ikdS  -6^ 

(Sk)c 

*  If  nT[(*(1>v)  V4>jl)  -n  *  (v<»-V«f)  -N(1)]dS  -6^ 

s. 

*  //NT  f  dS 

s,  « 

vn  =  normal  velocity 

x1  =  the  first  order  displacements 


(32) 


In  the  above  equation,  ne  and  Se  denote  that  corresponding  integrations  are  taken  over  a  single 

element  only.  The  velocity  potential,  at  boundary  nodes  of  a  single  element  consists  of 
three  components.  These  are  incident,  scattered,  and  radiated  potentials.  Hence,  the  element 
equations  can  be  written  as: 


[k]{4>?}  =  { f }j ;  j  =  H,S+D,i,k 


(33) 


Stiffness  Matrix  and  Nodal  Force  Vector  for  the  2D  Element 

Similarly,  the  element  equation  for  the  2D  case  can  be  derived  by  taking  the 
differentiation  on  the  functional  given  in  Equation  25  with  respect  to  •  Then: 
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ana) 


a*. 


a> 


=  {(VN)tF’(VN)  -4k2F*NTN}dxdy  -<{>| 


(2) 

e 


«2). 


-  /  F*  NTdS 


OV, 


dn 


•/4~’»*'(W(S) 


dS-4>f  =  0 


(34) 


and  a  system  of  equations  for  the  problem  is: 

[k(x,y>]  {❖f1}.  =  lUxj))^,  j=n,S+D,i,k 


(35) 


with 


[k(x,y)]  =  J J {(VN)tF*(VN)  -4k2F*NTN}dxdy 

-(sns)l“ 


(36) 


(f(x,y)l 


31 0 

d<PlI 

da 


NTdS 


(37) 


Equations  deducted  from  the  2D  and  3D  elements  can  be  assembled  to  form  a  set  of 
global  equations  pertaining  to  the  entire  fluid  region.  These  may  be  written  in  matrix  form  as: 

[K]  {4>a)}  =  |F)  <38> 


The  global  stiffness  matrix  [K]  is  symmetric  and  banded.  Its  half-band  width  indicates 
interactions  between  nodes  of  an  element  and  its  immediate  neighbors.  Equation  38  may  now 
be  solved  numerically. 
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SECOND-ORDER  WAVE  LOADS  ON  A  FLOATING  BODY 


Once  the  second-order  velocity  potential  is  determined,  the  second-order  wave  loading 
on  the  floating  body  can  be  calculated.  It  is  assumed  that  the  vessel  undergoes  small  first-order 

motions  X(1)  from  its  mean  position  X(0),  and  that  the  pressure  can  be  expanded  in  a  Taylor 
series  about  the  hydrostatic  pressure  at  the  mean  position  and  the  following  expression  is  found: 

p  =  p®  +  e  p®  +  e2  p®  +  O  (e3)  (39) 


where: 


hydrostatic  pressure:  p(0)  = 

-pgZ<°> 

(40) 

first-order  pressure:  p{,)  = 

-Pgz(1>  -  P*® 

(41) 

second-order  pressure:  p® 

=  -  i  p  |VO(1)|2  -  p  -  p(X0)  •  VO®) 

(42) 

The  total  hydrodynamic  force  exerted  on  the  body  relative  to  the  coordinate  system  Gxyz 
is  given  by  the  following  equation: 

F  =  -ffpNdS  (43) 

s 


where  S  is  the  instantaneous  wetted  surface  and  N  is  the  instantaneous  normal  vector  to  the 
surface  element  dS  relative  to  system  Gxyz . 

N  =  N®  +  eN(,) 


Substitution  of  the  pressure  p  given  by  Equation  39  and  retaining  the  integration  to  the 
second  order,  one  obtains: 


P  =  F(0)  +  ef0)  +  e2F0)  +  0(e3) 


(44) 


The  instantaneous  wetted  surface  S  is  composed  of  two  parts,  S0  and  s  (Figure  2).  S0  is  the 
constant  wetted  surface  up  to  the  still  water  line  and  s  is  the  surface  between  the  still  water  line 
and  the  wave  profile  along  the  vessel. 

The  second-order  forces  can  be  calculated  by  integrating  all  products  of  pressure  p  and 
normal  vector  which  give  second-order  force  contributions  over  the  wetted  surface  S0  and  s: 
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into  the  above  equation  yields  for  the  first  term 

-J/p®j¥®dS  =  6x-//p(l)fidS  =  0xF(I)  (46) 

So  **0 

where  F{1)  is  first-order  wave  exciting  force  and  hydrodynamic  reaction  force.  According  to 
Newton’s  second  law,  F(1)  can  be  written  as: 

Pa)  =  [M]^  <47) 

where  [M]  is  mass  of  vessel  in  the  air  and  Xq1}  is  the  first-order  acceleration  of  the  center  of 
gravity  of  the  body.  Thus,  the  first  term  of  the  second-order  force  may  be  obtained  from  the 
first-order  motion  dynamics  as: 

-ffpwN(VdS  =  0(1)  x  [Ml  X™  (48) 

So 

The  second  term  of  the  second-order  forces  is  integrated  numerically.  The  following  was 
obtained  from  Equation  45: 

-ff  p^fidS  =  p  //[-  |  V<t(,)t2  +  ®ta)  +  X(,)  •  V®!0  ndS  (49) 

sj2 
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The  integral  of  the  oscillating  surface  s  is  calculated  to  the  linearized  free  surface  by  substituting 
p(1^  from  Equation  41: 

-//p(1)SdS  =  -ff(~Pg  Z(1)  -  p  4>t(,))fidS  (50) 


where  the  surface  element  is  written  as: 

dS  =  dZ(1)  •  dl  (51) 


At  the  still  water  line,  the  linearized  free  surface  condition  to  first  order  is: 

-  p  >  =  p  g  C(1)  (52) 

where  is  the  first-order  free  surface  elevation.  Inserting  Equations  51  and  52  into  Equation 
50  results  in  the  following: 

-J/pO'HdS  =  -/•ipg(<ff#dl  (53) 

s  WL  2 

in  which  Cf  is  the  relative  wave  elevation  defined  by: 

=  C(1)  -  Z2}L  (54) 


The  final  expression  for  the  total  second-order  force  thus  becomes: 
p<2)  =  ~  f  j  Pg^ffiodl  +  6(l)  x  [M]  Xq> 


WL 


(55) 


SdS 
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The  second-order  force  given  in  Equation  49  can  be  divided  into  a  steady  component  and  a  time- 
dependent  biharmonic  wave  exciting  force.  The  steady  component  is  a  second-order  steady  wave 
drift  force  calculated  from  Equation  49  by  talcing  the  time  average  over  one  wave  period.  Since 
the  flow  is  periodic,  the  time  average  value  of  Cr  is  zero.  As  a  result,  the  second-order  steady 
wave  drift  force  becomes: 


-  /  Apg^ffidl  +  §(1)  x  [M] 

WL  2 


Jf  i  |v<&(1)|2  +  +  (x(1)  •  V4>t(1)) 

J  L  2 


ndS 


(56) 


For  regular  waves 


p  [ f  G^fidS  =  0 

So 


Thus,  the  second-order  force  can  be  written  as: 

FC)  =  ?D  ♦R.{f®e*1-1} 


(57) 


where  is  the  second-order  biharmonic  wave  force  and  is  given  by: 


?0)  «  ‘  /  “T"  (<1>1,))2  fi  dl  +  ®0>  x  [M]  xS) 
i  4g 


//[|p(VW  +  2iup4»(2)  +  2i o>  p(£(1) •  V4>(1)) 


(58) 


fidS 
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where: 


$0) 

=  d>«) 

e*"* 

=  4,(2) 

ei«t 

•OO) 

ag 

■  Xa 

e'"‘ 

e(1) 

=  0(1> 

e***' 

ca) 

n 

■e- 

<  x-S 

^  N* 

W 

ei«t 

The  total  hydrodynamic  moment  about  the  center  of  gravity  of  the  vessel  relative  to  the 
coordinate  system  Gxyz  is  given  by  Equation  5: 

M  =  - //p(X  x  N)dS  (59) 

s 


The  derivation  is  analogous  to  that  followed  for  the  forces.  The  final  expression  for  the  second- 
order  wave  moment  is: 

m®  =  -  /  |Pg (Cfi*  *  ♦  8®  *  me(1> 

WL  2 

(60) 


"//  '  “  p  l^*(1)l2  -  P  *?  ~  p(X(,)  * 


(X  x  n)  dS 


EXPERIMENTAL  STUDY 
Test  Setup  and  Procedures 

A  hydraulic  model  was  designed  to  observe  the  coupling  process  of  a  fixed  barge  with 
nonlinear  waves  in  shallow  water.  The  experiment  was  conducted  in  a  two-dimensional  wave 
tank  at  the  Civil  Engineering  Laboratory  of  Texas  A&M  University.  Figure  4  illustrates  the 
general  setup  of  the  experiment.  The  tank  was  36  meters  long,  0.91  meter  wide,  and  1.22 
meters  deep.  The  water  depth  was  maintained  at  a  constant  of  20.3  cm  throughout  the  entire 
experiment.  A  rectangular  box  (101.6  cm  by  33.9  cm  by  8.4  cm)  was  placed  at  roughly  11 
meters  from  the  wavemaker  and  centered  in  the  tank,  to  simulate  a  barge  afloat  at  the  water 
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surface.  The  box  was  mounted  to  a  rigid  aluminum  frame  with  three  legs  as  shown  in  Figure 
4.  The  box  was  further  ballasted  to  a  neutrally  buoyant  condition  at  the  draft  of  8.4  cm  in  the 
still  water 

Two  load  cells  were  placed  at  the  end  of  each  of  the  three  legs  to  measure  the  horizontal 
and  vertical  forces,  respectively.  Each  cell  provided  a  capacity  of  50  pounds  at  an  accuracy  of 
±1  percent  of  the  force  applied.  The  total  forces  and  moments  were  deduced  from  these 
measurements. 

A  significant  influence  on  the  wave  activity  around  the  test  model  due  to  the  close 
proximity  of  side  walls  of  a  narrow  wave  tank  is  anticipated.  The  wave  forces  imposed  on  the 
test  model  will  be  different  from  those  experienced  in  an  otherwise  open  water.  However,  this 
test  was  intended  for  validating  the  numerical  procedure.  The  side  walls  of  the  wave  tank  were 
simulated  in  the  numerical  procedure  to  give  a  fair  comparison  between  the  physical  and 
numerical  models. 

The  experimental  waves  were  generated  by  a  computer-controlled  Commercial  Hydraulics 
RSW  90-85  dryback,  hinged-flap  wavemaker.  The  water  pumped  by  the  flap  traveled  down  the 
tank  channel  and  evolved  into  shallow  water  Cnoidal  waves  before  it  reached  the  barge.  Two 
resistance-type  wave  gauges  were  used  to  monitor  the  wave  profiles.  One  was  placed  at  75  cm 
forward  of  the  bow  along  the  center  line  of  the  tank,  and  the  other  at  47  cm  aft  of  the  bow, 
centered  between  the  box  and  the  side  of  the  tank.  The  accuracy  of  the  gauges  was  ±  1  mm. 
Four  wave  periods  were  selected  intentionally  to  avoid  the  resonance  frequency  of  the  wave  tank. 
Six  wave  heights  in  the  range  that  promised  results  of  well  defined  Cnoidal  wave  profiles  were 
chosen  for  each  wave  period.  The  wave  parameters  used  in  the  test  are  summarized  in  Table  1. 
The  Ursell  number  of  these  waves  varied  from  0.5  to  70. 

Data  acquisition  was  accomplished  using  a  Hewlett  Packard  HP-3852A  data 
acquisition/control  unit  with  an  HP-330  workstation.  This  system  was  equipped  with  a  10- 
channel  relay  multiplexer  (HP-447 18 A)  and  an  integrating  voltmeter  operating  at  a  5-1 /2-digit 
resolution  (HP-44701 A).  Details  of  the  test  setup  and  procedures  are  described  in  Reference  14. 

Test  Results 

Twenty-four  sets  of  time  history  of  the  wave  profiles  and  the  horizontal  and  vertical 
forces  were  recorded.  A  sample  set  of  these  measurements  is  shown  in  Figure  5.  The  incident 
wave  elevation  shows  good  agreement  with  that  predicted  by  the  Cnoidal  wave  theory  in  both 
wave  height  and  period.  However,  it  is  noted  that  the  experimental  wave  is  asymmetric  in 
contrast  to  the  asymmetric  mathematical  form.  This  discrepancy  grows  as  the  wave  height 
increases.  This  results  from  the  wave  generation  mechanism  used  in  the  test.  The  force  time 
histories  are  similar  in  shape  to  that  of  the  wave  profile.  The  time  histories  of  all  measurements 
are  decomposed  by  Fast  Fourier  transform  up  to  the  third  harmonic  as  follows: 

3  X„ 

X(t)  =  X0  +  £  —  cos  (no>t  +  eB) 

B  -  1  2 
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where: 


Xq  =  the  mean  value 

Xjj  =  height  of  the  nth  harmonic  component 
en  =  the  phase  angle 
<i)  =  base  frequency 

The  symbols  H,  X,  and  Z  represent  the  wave  height  and  the  double  amplitudes  of  the 
forces  in  the  horizontal  and  vertical  directions,  respectively. 

It  was  found  that  all  the  mean  values  and  the  third  harmonics  were  negligibly  small  in 
comparison  to  the  first  and  the  second  harmonics,  and  are  therefore  not  discussed.  The  results 
of  decomposition  of  the  waves  are  shown  in  Figure  6.  The  height  ratios  of  each  harmonic 
component  to  the  resultant  wave  height  were  plotted  versus  the  resultant  wave  height.  The  same 
ratios  deduced  from  the  theoretical  Cnoidal  waves  are  presented  in  dashed  lines  in  the  same  plots 
for  comparison  with  the  measurements.  The  measurements  are  almost  identical  with  then- 
theoretical  counterparts.  It  is  noted  that,  for  the  short  period  waves,  the  height  ratios  remain 
essentially  constant  over  the  range  of  the  total  wave  height  considered.  For  the  long  period 
waves,  these  ratios  vary  nonlinearly  with  the  increase  of  total  wave  height  with  the  first  harmonic 
component  decreases  and  the  second  harmonic  increases.  This  result  clearly  demonstrates  the 
typical  nonlinear  characteristics  of  the  shallow  water  waves. 

The  force  time  histories  are  decomposed  likewise.  For  convenience,  the  components  are 
further  nondimensionalized  in  the  form  of  force  coefficients  as  follows: 

X. 

C  =  1 

*  pHjlb 

Z. 

c  *  1 

*  pHjlb 


where  p  is  the  specific  weight  of  water  and  "1"  and  "b"  are  the  length  and  the  beam  of  the  box, 
respectively.  These  force  coefficients  are  equivalent  to  the  ratios  of  the  heights  of  the  first 
harmonic  forces  to  the  first  harmonic  wave  height.  Figure  7  indicates  that  these  coefficients 
remain  constant  at  a  specific  wave  period.  This  implies  that  the  heights  of  the  first  harmonic 
forces  are  linearly  correlated  to  the  height  of  the  first  harmonic  wave.  The  second  harmonic 
terms  H2,  X2,  Z 2  are  plotted  as  functions  of  the  first  harmonic  wave  height  Hj  squared  as  shown 
in  Figure  8.  The  units  are  given  in  cm  for  H2,  Newtons  for  force  harmonics,  and  cm2  for  Hj2. 
The  results,  which  closely  fit  the  first-order  regression  line,  indicate  that  the  second  harmonic 
terms  X2,  Z2,  and  H2  are  linearly  correlated  to  the  square  of  the  first  harmonic  wave  height. 

The  proportionality  constants  for  X2  monotonically  increase  as  the  wave  period  increases, 
whereas  the  constants  for  Z2  generally  increase  but  slightly  fluctuate.  The  proportionality 
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constants  for  H2  increase  as  the  wave  period  increases,  although  the  constants  are  generally  very 
small  compared  to  those  of  the  forces. 

Nondimensional  Forces  and  Wave  Heights  Versus  the  UrseU  Number 

The  first  and  second  harmonics  of  the  forces  and  incident  wave  height  are  non- 
dimensionalized  by  dividing  them  by  their  respective  resultants,  then  compared  to  the  Ursell 
number, 


U. 


HL2 

h3 


where  h  =  water  depth.  The  first  harmonic  ratios  exhibit  a  linear  trend  and  the  second  harmonic 
ratios  present  some  nonlinearity,  as  shown  in  Figure  9.  Hj/H  and  Xj/X  decrease  linearly  as  the 
Ursell  number  increases  while  Zj/Z  remains  constant.  H2/H  and  X2/X  increase  initially  then 
approach  to  a  constant  value  while  Z2/Z  decreases  slightly  throughout  the  range  of  Ursell 
numbers.  This  trend  seems  to  agree  with  the  previous  finding  shown  in  Figure  10  from  Sarpkaya 
and  Isaacsson  (Ref  15)  on  nonlinear  shallow  water  wave  forces  on  a  vertical  circular  cylinder. 


NUMERICAL  RESULTS 

The  finite  element  procedure  derived  in  the  previous  section  on  theoretical  considerations 
was  executed  in  a  computer  program,  NAUTILUS,  coded  in  standard  Fortran  77  language.  This 
procedure  discretized  the  water  body  under  consideration  with  a  combination  of  three-  and  two- 
dimensional  finite  element  meshes.  This  feature  permitted  optimization  of  computer  resources 
while  retaining  the  essential  details  of  the  physical  problem  under  study.  The  discretization 
resulted  in  a  linear  matrix.  Equation  38,  which  was  treated  by  a  proven  procedure  based  on  the 
frontal  method  developed  by  Irons  (Ref  16).  The  solution  gave  a  complete  description  of  the 
fluid  activity  in  terms  of  velocity  potentials.  NAUTILUS  presents  the  results  in  a  convenient 
form  for  practical  application,  including  wave  height  distribution  over  the  entire  water  surface, 
wave-induced  forces  on  the  fixed  structures,  and  the  dynamic  motion  of  the  floating  structures. 
All  computations  are  executed  to  the  second  harmonic. 

NAUTILUS  was  used  to  reproduce  the  hydraulic  model  test  described  in  the  previous 
section  on  experimental  study.  The  results  were  used  to  verify  the  feasibility  of  the  NAUTILUS 
program.  Since  NAUTILUS  does  not  include  a  wave  generation  mechanism  to  simulate  the 
wavemaker,  the  wave  tank  will  be  approximated  with  a  channel  of  infinite  length  with  a  steady 
second-order  Stokes’  wave  approaching  from  one  end  of  the  channel.  The  water  body  within 
three  times  of  the  barge  length  from  the  center  of  the  barge  is  simulated  with  three-dimensional 
element  meshes,  while  the  rest  is  approximated  with  two-dimensional  element  meshes.  The 
water  body  under  consideration  extends  to  a  distance  of  ten  barge  lengths  on  each  side. 
Figure  1 1  illustrates  the  discretization  of  the  model  test  setup.  Wave  parameters  used  in  the 
model  test  are  entered  into  the  numerical  model  for  comparison  with  the  forces  observed  from 
the  model  test. 
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The  results  of  the  wave  forces  predicted  by  the  numerical  model  (in  black  symbols)  are 
presented  along  with  the  hydraulic  model  measurements  (in  white  symbols)  ir.  Figure  9.  The  fair 
agreement  between  the  theoretical  and  the  empirical  results  strongly  supports  the  approach  taken. 


CONCLUSIONS 

A  second-order  simulation  model  for  three-dimensional  wave-structural  couplings  in 
moderate  seas  has  been  formulated.  The  results  confirm  the  feasibility  of  reducing  a  weakly 
nonlinear  free  surface  problem  by  the  perturbation  method.  This  treatment  eliminates  the  need 
for  specifying  the  unknown  free  surface,  and  allows  the  nonlinear  problem  to  be  addressed  at  the 
well-defined  mean  water  surface  in  the  perturbation  components.  Each  component  forms  a  well- 
defined  linear  boundary  value  problem  that  can  be  solved  effectively.  The  present  study  solves 
the  component  boundary  value  problems  with  a  finite  element  procedure.  This  procedure,  which 
uses  a  combination  of  two-  and  three-dimensional  finite  elements,  gives  a  good  balance  between 
computation  efficiency  and  fair  proximation  of  the  physical  problem. 

The  wave-induced  forces  on  a  rectangular  box  held  steady  at  the  water  surface  were 
studied  both  numerically  and  experimentally  to  examine  the  significance  of  the  second-order 
forces.  The  second-order  forces  were  found  to  increase  with  the  Ursell  number,  LH2/h3,  which 
is  a  measure  of  the  wave  steepness  and  the  relative  water  depth,  in  a  manner  closely  resembling 
the  second-order  profile  of  a  Cnoidal  eave.  This  implies  that  the  significance  of  the  second-order 
forces  is  similar  to  that  of  the  second-order  wave  amplitude  to  the  total  amplitude  of  Cnoidal 
waves.  Hydraulic  model  measurement  indicates  that  the  second-order  lateral  forces  induced  by 
steep  waves  in  shallow  water  can  be  as  large  as  30  percent  of  the  total  lateral  force.  This 
confirms  the  importance  of  the  second-order  forces  in  high  sea  states.  Actual  wave  forces 
experienced  by  ocean  structures  in  rough  seas  can  be  much  greater  than  those  predicted  by  linear 
theories.  More  importantly,  the  additional  wave  forces  may  occur  at  frequencies  away  from  the 
design  wave  frequency,  and  coincide  with  the  resonance  frequency  of  the  structure  of  its 
components  and  cause  adverse  effects.  Ocean  structures  for  use  in  high  sea  states  must  be 
designed  to  withstand  second-order  forces. 
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Table  1 

Wave  Parameters  Used  for  the  Hydraulic  Model  Test 


T 

1  m 

(sec) 

(cm) 

0.75 

0.54,  1.09,  1.66,  2.11,  2.58,  3.13 

1.00 

0.52,  1.07,  1.60,  2.08,  2.55,  3.05 

1.20 

0.54,  1.52,  3.19,  5.09,  7.23,  9.23 

1.58 

0.69,  3.90,  4.98,  6.30,  8.48,  10.82 
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Figure  1 

Definition  sketch  of  a  compound  structure  and  the 
ambient  fluid  domains. 


2' 


Figure  5 

Profiles  of  sample  incident  wave  height,  horizontal  force, 
and  vertical  force. 
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Ratio  of  Wave  Height  Ra,io  of  Wavc  1  !e'ghl 
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Figure  7 

Hydraulic  model  measurements  of  the  wave  force  coefficients. 
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Order  Terms  (Newtons',  cm)  Second  Order  Terms  (Newtons,  cm ) 
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Figure  10 

Nonlinear  shallow  wave  forces  on  a  vertical  circular  cylinder. 
Ratio  of  Cnoidal  to  shallow  sinusoidal  wave  theory  force 
predictions  as  a  function  of  Ursell  number  U  for  various 
values  of  ka  (from  Ref  15). 


33 


Appendix 


EXTENSION  OF  BERKHOFF’S  WAVE  THEORY 
TO  THE  SECOND-ORDER  WAVE  PROBLEMS 


BASIC  EQUATIONS 


Field  equation: 


Free  surface  equation: 


where: 


v2<t>®  =  o 


1®1  a®  +  ii!? 

g  dz 


Q 

g 


=  z  =0 


Q  =  -i2w 


2  l  g 


Bottom  Condition: 


+  .  0> 

dx  dx  dy  dy  dz 


(A-l) 


(A-2) 


(A-3) 


The  derivation  of  the  reduced  equation  starts  with  basic  Equations  A-l  through  A-3  by 
introducing  dimensionless  coordinates  with  the  aid  of  the  characteristic  length  L  =  g/w2  at  the 
free  surface  and  the  mean  water  depth  H. 


A-l 


z 


L  =  wave  length 


H  =  mean  water  depth 
h(x,y)  =  water  depth 


Apply 


x 


/ 


_x 

L 


x  = 


X 

X 


y 


y  = 


y 

x 


and  drop  the  superscript  (2).  The  nondimensional  governing  equations  given  in  Equations  A-l 
through  A-3  are  given  as: 


a*<t>  ,  ,  ay  _  0 

dx'2  dy'2  \H/  dzn 


(A-4) 


ii  -  iigjjl  4,  =  f(<j>W)  =  ^  H 
dz'  g  g 


(A-5) 


ii  + 

H2 

ii  ii! 

+ 

ii!] 

=  0, 

z  =  -  h7 

(A-6) 

dz1 

XL 

dx'  dx 

+  dy' 

Qy  t 

A-2 


Introducing: 


H  .  o2H  H 

P  =  —  t  o  =  - ,  e  = 

L  g 


y/TL 


Equations  A-4  through  A-6  can  be  written  as: 

+  m  o  inQ 

dx'2  dyn  \i2  dz'2 


-  46<| >  =  f* ,  z'  =  0 
dz' 


i*  +  e2 

dz' 


(e±M+e±M]  m  0>  ,/..h/H..h/ 

, dx '  dx  dy'  dy  j 


Assuming  the  potential  function  4>  can  be  written  in  the  form: 

4>(x,y,z)  =  Z*(h/,pz,)<p(x/,y\ez') 

and  developing  the  function  <p  into  a  Taylor-series  with  respect  to  ez': 

4>(x,y>z)  =  Z*(<p0  +  c  z/ «Pj  +  e2z/2  <p2  +  •••) 


Assuming  a  mild  slope: 


z  = 


»  s  -  ^asp;- fs  «  i 

XL  ^  L 


Therefore: 


4>(x,y,z)  *  ZQi/>|iz/)*0(x/.y/) 


A-3 


(A-7) 


(A-8) 


(A-9) 


(A- 10) 
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Substitute  Equation  A- 10  into  A-7,  and  assuming  a  mild  slope  (e  <<  1) ,  we  obtain  a  simplified 
equation  as 


Z 


k  3x/2 


3z/2 


0 


That  is, 


<Po 


^<Po 

,ax/2 


i  a*z 

p2Z  3z'2 


The  left-hand  side  of  the  above  is  a  function  of  (x\y’)  only,  while  the  right-hand  side  is  a 
function  of  z’  only.  Therefore: 


<Po 


tax'2  dy' 


=  -4  k2 


(A- 11) 


and 

-i-JS.  .  4k*,  or  **  .  „*Z-4k* 
|i*Z  dzn  dz12 


Substituting  Equation  A- 10  into  A-9: 

Z  ♦  e*(  )  -  0,  z'.-h' 

dz 


Therefore: 


^<P0 

az; 


0  ,  z'=  -h',  e2-0 


(A-12) 


(A- 13) 


A-4 


Substituting  Equation  A- 10  into  A-8: 

—  <p-46Z<p=f*  or  —  -46Z  =  f*/<P0,  z'  =  0 
dz'  °  °  dz' 


Therefore,  the  simplified  equations  (assuming  e  <<  1)  are  given  as: 


Laplace  Equations: 


(&$  a2^ 


(dx/2  By'1 
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Free  Surface: 


Bottom  Condition: 


az 

dz' 


-  4  6  Z  =  f79c,  z'  =  0 


-^  =  0,  z'  —  -  h' 
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(A- 16) 


(A- 17) 


Solution  of  Z 


From  Equations  A-15  and  A-17,  the  following  result  can  be  obtained: 


Z(h',i')  =  cosh2k^(z'  ^h'> 
sinh4(k  ph') 


(A- 18) 
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or 


Z(h,z) 


cosh2f—  (z  +  h) 
\  L , 


sinh4(t h) 


Set  k/L  =  k  ...  wave  number 

Z(h,z) 


cosh2K(z  +  h) 
sinh4(ich) 


Solution  of  <p0 


From  Laplace  Equations  A-7  and  A- 10,  we  obtained: 
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From  Equation  A- 19,  the  integration  of  Z2  can  be  written  as: 
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The  group  velocity,  CB,  is: 
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From  free  surface,  bottom  condition  Equation  A- 16, 
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Therefore,  Laplace  Equation  A-20  becomes: 
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Equation  A-23  is  the  reduced  diffraction-refraction  equation. 
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